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ABSTRACT 


J 


A  linear  model  of  the  vertical  dispersion  of  near- 


inertial  waves  is  developed.  A  porosity  distribution  near 
the  bottom  of  the  computational  domain  minimizes  bottom 
reflections,  and  simulates  an  ocean  of  infinite  depth.  The 
model  is  used  to  show  that  the  vertical  dispersion  of  near- 
inertial  waves  in  the  upper  ocean  is  sufficient  to  account 
for  the  observed  rapid  decay  of  inertial  oscillations  in  the 
surface  layer.  When  the  upper  pycnocline  is  sufficiently 
peaked,  a  resonant  frequency  interference  effect  is 
predicted.  This  effect  modulates  the  dissipation  of  surface 
layer  inertial  oscillations,  and  their  magnitude  after  a 
storm  need  not  decay  monoton i cal ly.  We  look  at  deep-ocean 
observations  taken  during  MILE,  and  find  some  suggestive 
evidence  of  this  interference  effect.  We  also  compare  model 
simulations  with  Baltic  Sea  observations  made  by  KRADSS 
(1981),  and  show  that  near-inertial  waves  reflect  off  the 
shallow  (105  m)  bottom  within  a  few  inertial  periods  after  a 
storm.  ^ — 
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Section  1 
INTRODUCTION 


1 . 1  BACKGROUND 

Records  of  currents  in  the  upper  ocean  are  often 
dominated  by  intermittent  oscillations  at  the  inertial 
frequency.  It  is  well  known  that  the  strongest  of  these 
oscillations  are  wind  driven.  A  sudden  change  (over  a  time 
scale  less  than  f'1)  in  the  local  wind  stress  will  generate 
inertial  oscillations  in  the  ocean's  near-surface  layer. 
These  oscillations  then  decay  within  a  few  days.  The  dissi¬ 
pation  mechanism  is  not  yet  fully  understood.  In  this  report 
we  demonstrate  the  likelihood  that  vertical  dispersion  of  low 
frequency  internal  waves  is  responsible ,  at  least  in  part, 
for  the  dissipation  of  inertial  oscillations. 

Several  investigators  have  shown  that  inertial 
oscillations  in  the  surface  layer  are  largely  due  to  wind 
forcing.  POLLARD  and  MILLARD  (1970)  compared  the  inertial 
oscillations  computed  using  a  simple  horizontally  homogeneous 
wind-driven  model  with  upper  ocean  current  observations.  The 
action  of  the  wind  on  the  mixed  layer  was  imposed  as  a  body 
force.  To  dissipate  the  inertial  oscillations  and  to  obtain 
approximate  agreement,  it  was  also  necessary  to  introduce  an 
artificial  linear  damping  term.  Best  comparisons  were 
obtained  when  the  e-folding  decay  time  of  the  damping  term 
was  in  the  range  two  to  eight  days.  The  implication  is  that 
this  damping  term  is  a  gross  parameterization  of  a  certain 
process  (or  processes)  which  can  remove  energy  from  the  mixed 
layer. 


It  is  not  yet  clear  what  type  of  process  dominates 
the  dissipation  of  inertial  energy  in  the  surface  layer. 
During  a  wind  event,  strong  shear  near  the  inertial  fre¬ 
quency,  just  below  the  base  of  the  mixed  layer,  causes  the 
thermocline  interface  to  erode.  KRAUSS  (1981)  showed,  from 
observations  in  the  Baltic  Sea  and  from  numerical  results, 
that  the  Richardson  number  over  a  10  m  interval  across  the 
interface  is  driven  below  1/4  during  storms.  Mixing  due  to 
breaking  internal  waves  or  shear  instability  induces  entrain¬ 
ment  of  momentum  across  the  thermocline  interface.  The 
deepening  of  the  mixed  layer  distributes  momentum  over  a 
greater  volume,  and  inhibits  the  unlimited  growth  of  inertial 
oscillations. 

Nonlinear  interactions  may  help  to  dissipate 
inertial  waves.  KRAUSS  (1981)  cites  observations  of  inertial 
waves  below  the  Ekman  layer  with  amplitude  ~40  cm  s’*1,  and 
horizontal  wavelength  ~100  km.  Thus,  at  the  latitude  56*N, 
the  nonlinear  advection  term  uSu/8x  is  about  10%  of  the 
Coriolis  term  fu.  To  this  extent  then,  the  inertial  waves 
are  themselves  nonlinear.  The  resultant  super-inertial 
frequency  motions  generate  internal  waves  which  are  free  to 
propagate  downward  through  the  thermocline.  Spectral 
analysis  of  moored  current  meter  records  yields  some  evidence 
of  nonlinear  interactions.  Harmonics  at  twice  the  inertial 
frequency  (KRAUSS,  1981)  and  at  multiples  of  the  inertial  and 
semidiurnal  frequencies  (DAVIS,  DESZOEKE,  HALPERN,  and 
NIILER,  1981)  have  been  observed.  Further  studies  are  needed 
to  determine  the  effectiveness  of  these  nonlinear  inter¬ 
actions  in  removing  inertial  energy  from  the  surface  layer. 


Another  mechanism  for  removing  inertial  energy  was 
proposed  by  BELL  (  1978).  Undulations  in  the  base  of  the 
mixed  layer  associated  with  turbulent  eddies  are  advected  by 
horizontally  homogeneous  inertial  oscillations.  These 
motions  disturb  the  underlying  thermocline,  and  generate 
internal  waves  through  nonlinear  interactions.  The  transport 
of  momentum  by  the  radiating  waves  results  in  a  drag  force  on 
the  inertial  currents.  This  mechanism  may  contribute  to  the 
dissipation  of  surface  layer  inertial  energy,  but  so  far 
there  is  no  experimental  evidence  for  it. 

In  the  present  study  we  develop  a  linear  model  of 
vertical  dispersion  of  near-inertial  internal  waves.  The 
basic  mechanism  is  that  a  surface  mixing  layer  is  driven  by  a 
time  dependent,  spatially  varying  wind  stress.  Convergence 
of  near-inertial  currents  leads  to  Ekman  suction  (vertical 
velocity  field)  at  the  base  of  the  surface  layer.  The  Ekman 
suction  then  generates  internal  waves  which  propagate  down¬ 
ward  into  a  stratified  interior.  The  vertical  energy  flux 
due  to  these  waves  acts  as  a  radiational  damping  mechanism  on 
the  near-inertial  surface  currents. 

Energy  transport  is  proportional  to,  and  in  the 
same  direction  as,  the  group  velocity.  For  near-inertial 
internal  waves,  the  group  velocity  vector  is  directed  at  a 
very  shallow  angle  with  respect  to  a  horizontal  plane.  The 
vertical  component  of  group  velocity  increases  with 
increasing  stratification  and  horizontal  wavenumber.  The 
objective  of  this  study  is  to  determine,  for  realistic 
environmental  conditions,  whether  the  stratification  is 
sufficiently  strong  and  the  horizontal  wavelength  of  an 
imposed  wind  field  is  sufficiently  short  to  explain  the  rapid 
decay  of  inertial  oscillations. 


The  plan  for  this  report  is  as  follows.  In  Section 
1.2  we  present  the  basic  model  equations  of  motion.  In 
Section  2  we  perform  a  Fourier  series  expansion  in  the  x- 
direction,  and  describe  the  method  of  solution.  Section  3 
presents  a  validation  of  the  model  against  an  analytic 
solution  for  the  simple  special  case  of  an  unstratified 
diffusive  fluid.  The  dispersion  of  near  inertial  waves  into 
a  uniformly  stratified  fluid  is  described  in  Section  4.  The 
influence  of  stratification  on  the  dissipation  of  surface 
layer  inertial  oscillations  is  also  discussed.  In  Section  5 
we  consider  the  case  of  a  more  realistic  nonuniformly 
stratified  fluid.  We  compare  these  results  with  shallow- 
water  current  meter  observations.  In  Section  6  we  summarize 
the  results  of  our  work  during  FY81  and  FY82,  and  outline  our 
recommendations  for  future  work. 

1 . 2  EQUATIONS 

KROLL  (1975)  developed  a  wave  packet  propagation 
model,  in  which  the  Ekman  suction  at  the  base  of  a  viscous 
boundary  layer  drives  an  inviscid  stratified  interior.  The 
boundary  layer  is  decoupled  from  the  interior,  and  the  super- 
inertial  frequency  differential  of  the  interior  forcing  (the 
difference  between  the  Ekman  suction  frequency  and  the  local 
inertial  frequency)  is  prescribed  to  be  constant.  For  our 
purposes,  this  model  must  be  extended.  We  do  not  decouple 
the  boundary  layer  from  the  interior,  because  we  wish  to 
study  the  mutual  interaction  between  the  two  layers.  Also, 
we  do  not  specify  a  single  frequency  component,  but  treat  all 
frequency  components  of  the  forcing  simultaneously. 

In  certain  respects,  the  model  we  develop  is 
similar  to  that  of  KRAUSS  (1978,  1979,  1981).  Krauss'  model 


is  based  on  the  linearized  hydrodynamical  equations  of  motion 
in  a  viscous,  Boussinesq,  rotating  fluid.  Internal  gravity 
waves  are  generated  in  a  flat  bottomed  channel  of  finite 
width  and  infinite  length.  The  major  change  we  make  to 
KRAUSS'  model  is  the  addition  of  porous  damping  in  a  thin 
layer  just  above  the  bottom  boundary  of  the  computational 
domain.  The  purpose  of  this  porous  region  is  to  absorb 
internal  waves,  and  prevent  them  from  reflecting  off  the 
bottom  boundary.  We  also  extend  the  interpretation  of  the 
system  to  the  open  ocean. 


We  also  make  some  other  minor  changes  to  KRAUSS* 
model  equations.  We  ignore  the  horizontal  eddy  diffusion 
terms,  and  since  we  are  concerned  with  near- inertial  motions, 
we  assume  the  hydrostatic  approximation.  We  choose  an  ortho¬ 
gonal  coordinate  system  (x,y,z),  with  z  directed  positive 
upwards.  Table  1.1  is  a  list  of  the  variables.  For  a 
complete  development  of  the  model  equations,  we  refer  the 
reader  to  RUBENSTEIN  (1980).  The  equations  are  as  follows: 
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Table  1 . 1 

DEFINITIONS  OF  VARIABLES 


Right-handed  coordinate  system ,  z  positive  up, 
z  >  0  at  bottom,  D  at  surface. 

Velocity  components 

Time 

P'/Po*  where  p'  is  pressure  fluctuation  from 
a  reference  state,  and  pQ  is  a  representative 
value  of  density. 

Eddy  diffusivity,  a  function  of  z  only 
vaisala  frequency; 


pr(z)  is  a  reference  state  of  density. 

Buoyancy;  -p'g/p©,  where  p'  is  density  fluc¬ 
tuation  from  a  reference  state. 

Inertial  frequency  ■  2fisin(latitude) 

Porous  damping  coefficient;  zero  except  in  a 
thin  region  near  the  bottom  boundary. 


Here,  the  terms  ou  and  ov  are  the  damping  terms,  which  we 
discuss  further  in  Section  2.1.  The  N(z)  and  u(z)  distribu¬ 
tions  are  constant  in  time,  horizontally  homogeneous,  and 
represent  statistical  averages.  In  Section  4,  we  specify 
that  they  not  overlap  one  another.  If  there  were  significant 
turbulence  just  below  the  mixed  layer,  the  seasonal  thermo- 
cline  would  erode  further  and  entrain  additional  fluid  into 
the  mixed  layer.  Therefore  we  have  neglected  the  buoyancy 
eddy  diffusivity  term,  and  the  processes  which  maintain  the 
mean  buoyancy  profile  against  diffusion.  The  boundary 
conditions  are  as  follows: 


w  ■  0  at  z  «  0,  D, 


(1.6) 


“  0  at  z  “  °' 


•*  (!i  •  H)  •  •  at  2 
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The  surface  stress  vector  _r(x,y,t)  is  a  pre¬ 
scribed  function  which  drives  the  system.  Equation  (1.6) 
states  that  velocities  normal  to  the  top  and  bottom  boundary 
surfaces  should  vanish.  Equation  (1.7)  prescribes  that  the 
stress,  and  therefore  the  shear  at  the  bottom  should  be  zero. 
This  boundary  condition  allows  slippage  at  the  bottom  of  the 
computational  domain  which,  after  all,  is  an  artifice  of  the 
model . 


Section  2 
SOLUTION  APPROACH 


2.1  TRANSFORMED  EQUATIONS 

In  this  section  we  present  the  model  equations,  and 
develop  their  method  of  solution.  We  assume  that  the  wind- 
induced  surface  stress  vector  can  be  described  as  a  Fourier 
transform 


T(x,yft)  «  J  dk  exp(  ik  *x)j(k,t)  . 


(2.1) 


We  consider  the  wind-induced  stress  associated  with 
a  propagating  frontal  system.  We  choose  the  coordinate 
system  so  that  the  x-coordinate  is  aligned  in  the  direction 
of  propagation  of  the  front.  We  assume  that  the  front  is 
invariant  in  the  y-direction,  normal  to  the  direction  of 
propagation. 


As  our  model  equations  are  linear,  the  system 
separates  for  each  wavenumber  vector  k.  With  the  above 
assumptions  in  mind,  each  vector  k  is  directed  in  the  x- 
direction,  the  direction  of  front  propagation.  Writing 
k  ■  I  k I ,  the  general  solution  is 


/  dk  exp(ikx) 
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Dropping  the  hats,  the  equations  of  motion  become 
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and  the  boundary  conditions  are 


w  *  0  at  z  *  0,  D 
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(2.4a-c) 


The  shape  and  magnitude  of  o(z)  is  crucial  to  the 
reduction  of  bottom  reflections.  The  magnitude  of  o(z)  must 
be  sufficient  to  absorb  wave  packets  as  they  propagate  down¬ 
ward  through  the  porosity  region,  reflect  off  the  bottom,  and 
propagate  upward  again.  If  the  vertical  gradient  of  a(z)  is 
too  big,  then  waves  will  reflect  off  the  porosity  distri¬ 
bution  itself,  before  they  are  sufficiently  damped.  We  have 
a  great  deal  of  experience  using  porous  damping  formulations 
with  another  internal  waves  code,  GORWAK  (ROBERTS  and 
RUBENSTEIN,  1981).  A  useful  formulation  for  o(z)  is 
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o(z)  *  A  exp  [-sinh(sinh  z/B)  ] 


(2.5) 


The  constant  B  is  a  depth  scale ,  and  characterizes  the 
thickness  of  the  porosity  distribution.  The  gradient  do/dz 
decreases  with  increasing  B,  so  by  making  B  large ,  we 
minimize  reflections  off  the  porosity  distribution.  As  we 
increase  B  we  decrease  the  size  of  the  porosity-free  compu¬ 
tational  domain ,  so  we  strike  a  compromise,  and  let  B  *  0.4D. 
Through  experimentation,  we  find  that  A  ~  f  is  an  optimal 
value  which  minimizes  reflections  off  the  bottom. 

2.2  METHOD  OF  SOLUTION 


2.2.1 


Finite-Difference  Mesh 


We  choose  to  solve  the  equations  and  boundary 
conditions  (2.3),  (2.4)  using  a  finite  difference  scheme. 
The  scheme  is  leap  frog  in  time,  and  staggered  and  stretched 
in  depth.  Figure  2.1  shows  a  schematic  representation  of  the 
finite  difference  mesh.  The  time-step  constraint  due  to 
numerical  instability  of  the  diffusion  terms,  namely 


6t  <  (5z)2/n, 


(2.6) 


is  very  serious.  With  6z  «  100  cm  resolution,  and  a  typical 
value  of  eddy  diffusivity  m  ■  100  cm2/sec,  the  time  step  6t 
must  be  less  than  100  sec.  In  order  to  avoid  the  constraint 
(2.6),  we  construct  a  hybrid  implicit  scheme,  where  the  dif¬ 
fusion  and  porosity  terms  are  implicitly  formulated.  The 
other  terms  remain  explicit. 
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Figure  2. 
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1  Schematic  of  finite  difference  mesh.  The  mesh  is 
staggered  in  depth  and  in  time.  The  mesh  is 
divided  into  n  unequal  intervals  in  depth.  The 
variables  enclosed  in  parentheses  are  outside  the 
domain ,  and  are  determined  by  the  boundary 
conditions. 


The  depth  coordinate  is  stretched  according  to  the 


formula 


In  (  y/(  1- 


.in 


(2.7) 


where  j  is  an  integer  or  half-integer  index,  n  is  the  number 
of  mesh  intervals,  and  y  is  a  stretching  parameter.  We  can 
solve  for  z  to  get 

z/D  -  1  +  Y  -  exp  [iji  In  ( y/(  1  +  Y)  )  +  ln(l+Y)  ].(2.8) 

2.2.2  Time-Step  Increment 

The  predictive  equations  (2.3a),  (2.3b),  and  (2.3d) 
for  u,  v,  and  b  are  integrated  by  computing  the  increments 
6u,  6v,  and  6b,  and  accumulating  these  increments  to  advance 
forward  to  the  next  time  step.  For  example,  we  introduce  our 
notation  for  the  v-momentum  equation; 

vi+1  -  v1  +  6v ,  (2.9) 

where  the  superscripts  indicate  the  time  step.  We  write  the 
v-momentum  equation  in  the  form 

vi+1-  v1  -  -ful  +  5_iiS-(vi  +  96v)  -  c(vl  +  66v),  (2.10) 
- 5t -  1  z 

where  0  is  a  parameter  which  specifies  the  degree  of  implic¬ 
itness.  For  all  cases  discussed  in  this  report,  we  let 
0  *  1/2.  Rearranging,  (2.10)  becomes 

(--1  +  ©o  -  06  ud..)6v  *  -fu1  +  a  uSv1  -  ov*.  (2.11) 


The  expression  in  parentheses  on  the  left  hand  side  may  be 
represented  by  a  tridiagonal  matrix.  We  solve  for  &v  using  a 
vectorized  tridiagonal  solving  routine  written  by  J.  Boris 
for  use  on  the  Texas  Instruments  Advanced  Scientific  Computer 
at  the  Naval  Research  Laboratory. 

2.2.3  Momentum  Transport  Constraint 

We  can  integrate  the  hydrostatic  balance  equation 
(2.3c)  to  obtain  the  pressure  to  within  a  constant,  a; 

P  *  Iq  bdz  +  a  .  (2.12) 

We  split  the  time-step  increment  6u  into  two  components, 

6u  *  6u  +  a6u,  (2.13) 

where  6u  and  a6u  are  driven  by  the  two  components  of  (2.12), 
respectively. 

The  first  component  6lT  is  the  solution  to  the 
finite  difference  equation 

(-5^  +  So  -  ®&24&z)  6u  -  fvA  -  /*  bdz 

+  a  ua  u1  -  ou1,  (2.14) 

Z  Z 

where  the  buoyancy  integral  gives  the  partial  pressure,  apart 
from  a.  Also,  6TT  satisfies  the  inhomogeneous  boundary  condi¬ 
tion  at  the  surface 


z 


and  the  homogeneous  boundary  condition  at  the  bottom 


^8  (  6u)  «  0;  z  *  0. 
z 


(2.16) 


The  second  component  6u  is  the  solution  to  the 


equation 


fat  +  60  “  *V&z  6“  * 


(2.17) 


and  satifies  the  homogeneous  boundary  conditions 


Ud_( 6u)  «  0;  z  ■  0 f D. 
z 


(2.18) 


Note  that  the  expressions  in  parentheses  in  (2.11),  (2.14) 
and  (2.17)  are  identical.  Also  note  that  (2.17)  and  (2.18) 
are  independent  of  time,  so  we  need  to  solve  for  6u  only 
once.  In  order  to  solve  for  a  we  impose  an  additional 
constraint 


fn  udz  *  0  r 


(2.19) 


obtained  by  integrating  the  continuity  equation  (2.3e)  with 
respect  to  depth  and  applying  the  boundary  condition  (2.4a). 


2.2.4 


Solution  Algorithm 


In  this  section  we  outline  the  steps  of  the 
solution  algorithm. 


Initialize  uf  v,  b,  and  p  to  zero.  Multiply  N2 
by  the  chosen  k2  value. 


Using  (2.17)  and  (2.18),  solve  for  6u. 

Integrate  the  hydrostatic  relation  (2.3c)  to 
obtain  the  (partial)  pressure; 

z 

p(z)  -  /  b  dz.  (2.20) 

0 

Solve  (2.11)  for  6v ,  subject  to  the  boundary 
conditions  (2.4).  Use  (2.9)  to  get  the  new 
updated  value  of  v. 

Solve  (2.14)  -  (2.16)  for  6u. 

Perform  a  preliminary  update  on  u  to  get  an 
intermediate  value; 

u*  «  u  +  6u.  (2.21) 

Use  the  constraint  (2.19)  to  solve  for  a; 

D  * 

”/  u  dz 


/  6u  dz 
0 


Note  that  the  denominator  in  (2.22)  is  constant 


B.  Perform  an  update  on  u; 

u  *  u*  +  a6u.  (2.23) 

9.  Integrate  the  continuity  equation  (2.3e)  to 
obtain  w; 

z 

w(z)  »  /  u  dz.  (2.24) 

0 

10.  Perform  an  update  on  b; 


b  ■ 

new 

bold  ' 

6t  YT  N^w 

(2.25) 

Repeat 

steps 

(3)  -  (10) 

until  the  desired 

integration  time  is  reached. 


Section  3 
ODE  VALIDATION 


In  this  section  we  validate  the  numerical  code  for 
the  special  case  of  an  unstratified  fluid,  where  an  analytic 
solution  is  available.  The  equations  can  be  simplified  if  we 
-define  a  new  time  variable  t'  *  tf  and  a  new  eddy  diffusivity 
V  :  a.-uf1.  This  allows  us  to  eliminate  the  Coriolis 
parameter  f  from  the  equations.  We  set  N  ■  a  *  0  in  (2.3), 
and  drop  the  primes,  leaving  us  with  the  equations 

afcu  -  v  s  dz(  ^dzu>  * 

5fcv  +  u  -  dz(  ^azv)  * 

W  -  k2u  =  0  .  (3.1a-c) 

We  examine  the  special  case  where  the  eddy  diffusivity  is 
-constant r  ' p.  ■  10~3.  Note  that  since  stratification  has  been 
-set  to  zero,  the  buoyancy  and  pressure  fluctuations  hare  left 
-the  problem.  The  continuity  equation  (3.1c)  relates  w  and  u, 
but  -w  is  a  passive  variable,  so  it  can  be  safely  ignored.  We 
let  V  •  u  +  iv,  so  that  (3.1a,b)  can  be  combined  to  read 

(  ?.  +  i  ”  li*  „)U  *  0  .  (3.2) 

c  ZZ 

For  boundary  conditions  we  let 

ui  C  =  x  z  ■  0,  (3.3) 

u  -  C  z  ♦  -•  ,  (3.4) 

where  t  is  the  complex  surface  stress.  ; 


From  GONELLA  (1971),  the  general  solution  starting 
from  an  initial  velocity  UQ(z)  *  U(z,0)  is 

U(z,t)  -  U.(z,t)  *  u_(z)  +  U.(z,t)  *  2t(t),  (3.5) 


where 


U.(z,t)  *  Y(t)  - - - —  exp  ( -  z  2/4  |it )  , 

1  (4xut)1/2 


(3.6) 


and  Y(t)  is  the  Heaviside  step  function.  The  asterisks 
denote  convolution  operations,  the  first  with  respect  to  z 
and  the  second  with  respect  to  t. 


We  consider  the  special  case  where  the  initial 
velocity  is  zero,  and  the  x-component  of  surface  stress  is 
impulsive; 


UQ(z)  -  0  , 


(3.7) 


T(t)  «  6(t) 

Substituting  (3.6)-(3.8)  into  (3.5)  we  get  U(z,t) 
2U-|(z,t).  At  the  surface  z  *  0,  the  solution  is 


U(0,t) 


Y(t)  e 


(  «»it)  1/2 


(3.9) 


The  solution  describes  inertial  oscillations  with  an 
amplitude  which  decays  as  t”1/2.  The  reason  for  this  decay 
is  the  diffusion  of  momentum  through  a  surface  layer  whose 
thickness  grows  as  t1/2. 
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We  integrated  (3.1)  using  the  numerical  code  for  8 
inertial  periods.  The  number  of  mesh  intervals  was  n  *  90, 
and  the  stretching  parameter  in  (2.7)  was  y  -  0.1,  which 
corresponds  to  mesh  spacing  exactly  eleven  times  finer  at 
z  =  1  than  at  z  »  0.  The  time  step  6t  was  x/32.  In  Figure 
3.1,  the  first  2.5  inertial  periods  of  the  numerical  solution 
are  compared  with  the  analytic  solution  (3.9).  Because  the 
depth  average  of  u  in  the  numerical  model  is  constrained  to 
be  zero,  and  the  bottom  boundary  condition  allows  slip,  we 
subtracted  the  bottom  velocity  from  the  surface  velocity 
before  displaying  the  results.  The  numerical  solution  has 
been  corrected  for  the  lag  which  occurs  because  the  surface 
stress  is  not  a  true  impulse,  but  is  applied  over  the 
duration  of  the  first  time  step.  With  these  corrections,  the 
comparison  is  excellent. 

Figure  3.2  gives  the  reader  a  good  feel  for  the 
overall  behavior  of  the  solution.  In  this  figure,  the  v- 
component  of  velocity,  computed  numerically,  is  plotted  as  a 
function  of  z  and  t.  The  first  3.75  inertial  periods  are 
shown.  One  can  clearly  see  the  inertial  oscillations,  which 
are  strongest  in  the  surface  layer.  The  depth  of  the  surface 
layer  is  delineated  by  the  gently  sloped  solid  line  which 
represents  zero  amplitude.  The  depth  of  the  surface  layer 
increases  as  t1/2  for  small  t.  For  large  t,  the  depth 
asymptotically  approaches  z  *  1/2,  for  the  reasons  described 
above,  in  the  previous  paragraph. 
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Fig.  3.1  Comparison  between  analytic  solution  at 
surface  (solid;  u  component  and  dashed; 
v  component  curves)  and  numerical  solution 
(crosses  and  diamonds) .  An  impulsive 
wind  stress  is  applied  to  the  surface  of 
an  unstratified,  rotating  fluid  at  time 


Section  4 

INERTIAL  WAVES  IN  A  SURFACE  LAYER  OVERLYING 
A  UNIFORMLY  STRATIFIED  THERMOCLINE 


In  this  section  we  examine  the  dispersion  of  near- 
inertial  internal  waves  in  a  uniformly  stratified  thermo- 
cline,  and  their  resultant  decay  in  the  surface  layer.  We 
consider  here  only  impulsive  surface  stress  forcing.  This 
rather  extreme  simplification  of  real  wind  events  allows  us 
to  separate  out  the  effects  of  phase  interference  of  surface 
layer  inertial  oscillations ,  which  accompany  wind  events  of 
finite  duration. 


We  examine  the  idealized  case  where  both  strati¬ 
fication  and  eddy  diffusivity  are  step  functions  in  depth; 


N(z) 


N0  0  <  z  <  h 
0  h  <  z  <  D 


(4.1) 


!0  0  <_  z  <  h 

(4.2) 

MO  h  <  Z  <  D  . 

The  depth  D-h  is  the  mixed  layer  depth.  The  eddy  diffusivity 
Mo  allows  momentum  to  diffuse  throughout  the  unstratified 
region  h  £  z  <  D.  The  N(z)  and  ii(z)  distributions  do  not 
overlap.  Momentum  cannot  diffuse  below  z  <  h,  but  must  be 
transported  via  internal  waves  through  the  stratified 
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!? 
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interior.  Table  4.1  lists  the  parameters  for  all  of  the 
cases  discussed  in  Section  4. 

In  all  cases  in  this  section  and  in  Section  5,  the 
time  step  increment  was  1/64  of  an  inertial  period ,  or  982 
sec.  The  surface  stress  driving  function  was  an  impulse  in 
the  y-direction.  This  impulse  lasted  for  the  duration  of  the 
first  time  step,  and  was  then  turned  off.  The  reason  for 
applying  an  impulse  is  that  we  wanted  to  study  the  effects  of 
vertical  dispersion  in  isolation  from  other  factors.  If  we 
had  used  a  realistic  wind  event,  lasting  an  appreciable 
fraction  of  an  inertial  period,  then  phase  interference 
effects  would  become  important. 

For  a  similar  reason,  we.  chose  values  of  eddy 
diffusivity  which  are  probably  unrealistically  large.  The 
role  of  diffusivity  is  simply  to  redistribute  momentum 
through  the  surface  layer.  Tests  show  that  the  depth- 
averaged  energy  in  the  surface  layer  and  in  the  stratified 
region  is  quite  insensitive  to  the  value  of  nQ.  The 
advantage  to  using  a  large  diffusivity  is  that  the  momentum 
redistribution  can  occur  quickly.  Therefore,  the  dissipation 
of  surface  layer  inertial  oscillations  due  to  internal  wave 
dispersion  is  separated  from  that  due  to  diffusion. 

* 

The  magnitude  of  the  surface  stress,  applied  during 
its  initial  impulse,  was  0.08  m2s”2.  If  this  stress  were 
distributed  over  a  realistic  "frontal  passage  time",  say,  1/4 
of  an  inertial  period  (about  4.4  hours),  then  this  stress 
impulse  is  equivalent  to  a  distributed  magnitude  of  5.3x10"** 
m2s"2.  Therefore,  the  applied  surface  stress  is  equivalent 
to  that  due  to  a  sudden  wind  event  lasting  4.4  hours,  and 
having  a  magnitude  of  19  ms**1  (about  37  kt). 
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Table  4.1 

Parameters  for  Cases  in  Section  4 


Mixed  Layer  Depth 


NQk(  s^nT1 ) 

d(m) 

H 

i 

a 

C4 

o 

n 

Porosii 

1 .25x1 0-6 

20 

0.02 

Yes 

1 .25x1 0”6 

20 

0.02 

No 

1 .25x10“6 

80 

0.08 

Yes 

0 

20 

0.02 

NO 

3.1 25x 1 0“7 

20 

0.02 

Yes 

6.25X10'7 

20 

0.02 

Yes 

2.5x10“6 

20 

0.02 

Yes 
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EFFECTS  OF  BOTTOM  POROSITY 


The  principal  effect  of  bottom  porosity  is  to 
minimize  reflections  of  internal  waves  off  the  bottom 
boundary.  In  order  to  demonstrate  this  effect  we  performed 
two  model  simulations,  Case  1  with  porosity  A  *  1.5x10“** 
s”1  ,  B  *  320  m  in  equation  (2.5))  and  Case  2  without 
porosity.  The  values  of  N0k,  d  ■  D-h,  and  n0  in  (4.1) 
and  (4.2)  were  1.25x10”6  s^m”1  ,  20  m,  and  0.02  m2s“1, 
respectively.  The  computational  bottom  depth  D  was  800  m, 
the  horizontal  wavelength  was  100  km,  and  the  Vaisala  fre¬ 
quency  was  5.7  cph. 


After  integrating  for  16  inertial  periods,  we  com¬ 
puted  the  Fourier  transforms  of  the  v-component  of  surface 
velocity,  shown  in  Figure  4.1.  The  simulations  were  nearly 
identical  for  low  frequencies  «  <  f.  Both  simulations  show  a 
strong  resonant  spectral  peak  at  u  *  f.  Above  f,  the  porous 
simulation  spectrum  decays  monotonically ,  while  the  porosity- 
free  case  has  a  series  of  spectral  peaks.  We  can  understand 
these  peaks  if  we  consider  the  dimensional  dispersion 
relation  in  a  rotating  hydrostatic  fluid; 


k2N2 


(4.3) 


where  <*>  is  angular  frequency,  and  k  and  A  are  the  horizontal 
and  vertical  components  of  the  wavenumber  vector.  In  terms 
of  the  Vaisala  frequency  N0,  the  dispersion  relation  for 
the  j'th  vertical  mode  is 


0) 

M 

3 

•H 

En 


u>  . 

3 


(4.4) 


2  »  f2  +  ( NQkD/jn) 2  . 

The  normalized  frequencies  of  the  lowest  four  modes  for  Case 
2,  with  N0kD  *  10"3  s_1 ,  are  wj/f  ■  3.34,  1.88,  1.46, 
1.28;  j  *  1,  2,  3,  4.  These  frequencies  correspond  to  the 
four  spectral  peaks  in  Figure  4.1. 

Thus,  we  see  that  the  bottom  porosity  distribution 
effectively  eliminated  reflections  from  the  bottom  boundary, 
and  prevented  the  establishment  of  vertical  modes  of  internal 
waves.  RUBENSTEIN  (1981,  1982)  examined  the  porous-free  case 
in  detail. 


4.2  VERTICAL  DISPERSION 

In  this  subsection  we  examine  the  vertical  dis¬ 
persion  of  internal  waves  in  greater  detail.  We  include  a 
bottom  porosity  distribution  in  order  to  minimize  bottom 
reflections. 

Figure  4.2  shows  the  response  of  the  v-component  of 
velocity  to  a  surface  stress  function  applied  impulsively  at 
time  t  »  0  in  the  x-direction,  for  Case  3.  The  values  of 
N0k,  d,  and  in  (4.1)  and  (4.2)  are  1.25x10-6  s"1m”1, 

80  m,  and  0.08  m2s"1,  respectively.  The  displayed  z  domain 
is  unstretched,  and  ranges  from  the  surface  down  to  a  depth 
of  520  m.  Below  that  depth  the  solution  is  not  shown, 
because  of  the  porous  damping  layer  there.  The  time  axis 
range  from  t  ■  0  to  t  *  7.5*/f,  or  3.75  inertial  periods. 


-6 


Within  half  an  inertial  period  after  t  =  0, 

momentum  becomes  fairly  well  diffused  through  the  surface 
layer,  and  within  a  full  inertial  period  (t  =  2n/f),  the 
shear  in  the  surface  layer  approaches  zero.  Inertial  oscil¬ 
lations  are  very  strong  in  the  surface  layer,  and  the 
strength  slowly  decays  in  time.  Meanwhile,  internal  waves 
propagate  into  the  stratified  interior.  From  (4.3)  we  obtain 
the  vertical  component  of  group  velocity  cz  *  bu/bl; 


k2N  2 
U>A3 


(4.5) 


As  the  Fourier  transform  of  a  delta  function  is  a  constant, 
the  initial  wind  impulse  transmits  all  frequencies  equally  to 
the  fluid.  From  (4.3),  the  highest  frequencies  are  asso¬ 
ciated  with  the  smallest  A.  For  u  >>  f,  the  vertical  group 
velocity  is  approximately  c2  ~  kN/A2.  Thus  the  internal 
wave  components  with  the  smallest  A  and  the  largest  w  are 
able  to  propagate  downwards  most  rapidly,  away  from  the 
surface.  After  these  waves  propagate  away,  the  upper-thermo- 
cline  region  is  left  with  waves  with  near-inertial  fre¬ 
quencies  w  >  f.  These  waves  propagate  vertically  much  more 
slowly. 


In  the  upper  thermocline  the  vertical  phase 
velocity  w/A,  according  to  our  argument  above,  should 
decrease  with  time.  In  Figure  4.2  shaded  and  unshaded  areas 
correspond  to  regions  of  positive  and  negative  v-component  of 
velocity,  and  therefore  are  bands  of  constant  phase.  As  time 
progresses  these  bands  slope  progressively  closer  to  the 


4-8 


horizontal.  This  result  confirms  our  expectation  that 
vertical  phase  velocity  decreases  with  time,  and  so  too  does 
vertical  group  velocity. 

Another  feature  of  interest  is  the  velocity  maximum 
which  develops  at  the  top  of  the  stratified  interior  region. 
Strong  vertical  shears  develop  just  above  this  velocity  maxi¬ 
mum.  For  convenience,  let  us  say  that  at  some  instant  in 
time,  a  particular  vertical  wavelength  component  AQ  is 
dominant  in  the  upper  thermocline.  Longer  wavelength  compon¬ 
ents  A  >  A0  have  already  had  time  to  propagate  away. 
Shorter  components  A  £  A0  are  in  phase  at  the  density 
interface  and  constructively  interfere  just  below  the  density 
interface,  resulting  in  a  velocity  maximum  there.  As  depth 
below  the  interface  increases,  the  different  components  begin 
to  interfere  destructively.  At  some  later  time  the  A0 
component  will  have  propagated  away,  leaving  some  shorter 
dominant  component.  Hence  the  vertical  extent  of  the 
velocity  maximum  becomes  progressively  thinner  with  time. 
However,  the  amplitude  of  this  feature  does  not  significantly 
decay  over  the  course  of  simulation.  After  16  inertial 
periods,  the  velocity  maximum  decreases  in  strength  by  only 
25%. 

4.3  ENERGY  PARTITION 

We  have  examined  the  energy  balance  in  the  surface 
layer  and  in  the  stratified  region  immediately  below.  Figure 
4.3  shows  the  depth-averaged  and  time-averaged  (over  1/2 
cycle)  kinetic  energy  and  potential  energy  in  the  ranges 
0.975  £  z/D  £  1.0  (surface  layer)  and  0.95  £  z/D  £  0.975 
(upper  stratified  region)  for  Case  1.  Kinetic  energy  density 
is  given  by  p|u|2/2,  and  potential  energy  density  is 
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Time  (Inertial  Periods) 

3  Depth-Averaged  kinetic  energy  in  surface  layer 
(A) ,  in  upper  stratified  region  (B) ,  and 

potential  energy  in  upper  stratified  region 
(C)  ,  for  Case  1. 


pb2/2N2,  where  p  is  mass  density  and  b  is  buoyancy.  Since 
the  surface  layer  is  unstratified,  the  potential  energy  there 
is  zero.  Also,  since  the  vertical  length  scale  D  is  much 
smaller  than  the  horizontal  wavelength  \  «  2n/k,  the  contri¬ 
bution  of  vertical  velocity  to  the  kinetic  energy  is  negligi¬ 
ble. 


A  statement  of  energy  conservation  is  written 

dfc  [-j  p  ( j  u  |  2  +  b2/N2  )  ]  +  V*(pu>]  =  0.  (4.6) 


In  the  absence  of  energy  flux  pu,  kinetic  energy  can  only 
be  gained  (or  lost)  by  transforming  from  (or  into)  potential 
energy.  Since  a  fluid  particle's  buoyancy  magnitude  |b|  has 
two  maxima  and  minima  per  wave  cycle,  kinetic  energy  also  has 
two  maxima  and  minima  per  cycle.  For  this  reason  the  curves 
in  Figure  4.3  represent  averages  over  1/2-inertial  cycles. 
In  the  remainder  of  this  section,  we  represent  this  time- 
average  operation  by  a  set  of  angled  brackets. 


In  the  upper  stratified  region,  kinetic  energy  is 
greater  than  the  potential  energy.  From  LEBLOND  and  MYSAK 
(1978),  this  ratio  for  a  single-frequency  plane  wave  is 


<E_> 


<Ek> 


N  2  (  0)2-f  2  ) 

w2  (N  2-f  2  )  +  f  2  (N  2-  a)2  ) 


(4.7) 


The  brackets  denote  an  average  over  a  complete  cycle.  Since 
N  >>  a)  >>  f,  we  can  approximate; 


Note  that  this  ratio  is  independent  of  N,  and  becomes  zero  in 
the  limit  u  *  f.  Inverting  (4.8)  to  solve  for  w,  we  get 


«  *  f  1+<V/<V 


-<Ep>/<Ek> 


(4.9) 


If  we  keep  in  mind  that  our  numerical  solution  represents  a 
continuum  of  frequencies,  we  might  still  use  (4.9)  in  con¬ 
junction  with  the  results  in  Figure  4.3  to  estimate  the 
dominant  frequency  present  in  the  upper  stratified  region. 
Figure  4.4  shows  (w-f)/f,  where  u>  is  the  dominant  frequency 
in  the  depth  range  0.95  £  z/D  £  0.975.  We  see  that  the ‘fre¬ 
quency  differential  decays  approximately  as  t“3/2.  During 
the  time  period  1  £  t(f/2*)  £  5,  the  dominant  frequency  <*> 
decays  from  1.2f  to  1.015f.  This  result  is  consistent  with 
near-surface  observations.  Over  a  wide  variety  of  con¬ 
ditions,  GONELLA  (1971),  PERKINS  (1972)  and  KUNDU  (1976) 
observe  spectral  peaks  at  frequencies  ranging  from  1.007f  to 
1 . 1 1f . 


I 


It  is  possible  to  estimate  the  dominant  vertical 
component  of  the  group  velocity  in  the  upper  stratified 
region.  Combining  (4.3)  and  (4.5)  we  get 


y-f2V 

wkN 


(4.10) 


We  use  (4.9)  to  provide  an  estimate  of  the  dominant  fre¬ 
quency.  Figure  4.5  shows  cz  as  a  function  of  time,  in  the 
depth  range  0.95  £  z/D  £  0.975.  The  vertical  component  of 
group  velocity  cz  decays  approximately  as  t“2*4. 
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Figure  4.5  Vertical  component  of  group  velocity 
in  the  upper  stratified  region,  for 
Case  1. 


4.4 


STRATIFICATION 


One  of  our  primary  concerns  is  the  effect  of  strat¬ 
ification  on  the  dispersion  of  near-inertial  waves.  In  this 
section  we  examine  a  set  of  cases ,  each  with  a  different 
density  stratification,  and  determine  how  the  dissipation 
rate  of  inertial  oscillations  in  the  surface  layer  reacts  to 
varying  stratification. 

As  throughout  Section  4,  the  Vaisala  frequency  and 
eddy  diffusivity  profiles  are  of  the  forms  given  by  (4.1)  and 
(4.2).  We  compared  Cases  1,  4,  5,  6,  and  7,  described  in 

Table  4.1.  The  values  of  N0k  were  0,  3.125,  6.25,  12.5, 
and  25  (x10~6  s~1m~1).  The  mixed  layer  depth  d  was  20  m,  and 
Mo  was  0.02  m2s“1. 

The  purpose  of  Case  4,  with  N0k  *  0,  is  to 

describe  a  base  state  in  which  all  of  the  energy  is  confined 
to  the  surface  layer.  Since  the  unstratified  fluid  cannot 
sustain  internal  gravity  waves,  there  is  no  need  for  bottom 
porosity.  As  a  result,  momentum  is  conserved~the  net 
momentum  is  zero — but  energy  is  not  conserved.  Immediately 
after  the  initial  surface  stress  impulse,  a  sharp  momentum 
gradient  is  present  at  the  surface.  Diffusion  redistributes 
the  momentum  through  the  surface  layer,  and  quickly  makes  the 
velocity  distribution  uniform  with  depth  in  the  surface 
layer,  and  the  kinetic  energy  approaches  a  constant  value. 

In  Case  4,  the  kinetic  energy  approaches  within  It  of  its 

constant  final  value  after  time  t  ■  x/5f  (1/10  of  an  inertial 
period).  Therefore,  we  can  obtain  a  measure  of  the 
dissipation  rate  of  inertial  oscillations  due  to  vertical 
dispersion,  by  comparing  the  evolving  depth  averaged  kinetic 
energy  in  the  surface  layer  with  its  "post-diffusion"  value 
at  time  t  ■  x/5f. 


Figure  4.6  shows  the  time  evolution  of  depth 

averaged  kinetic  energy  in  the  surface  layer  for  Cases  1,  5, 
6,  and  7.  The  dissipation  rate  increases  with  increasing 

N0k.  For  comparison  purposes ,  the  e-folding  dissipation 

time  scale  is  a  useful  quantity,  and  is  tabulated  in  Table 

4.2.  These  time  scale  estimates  are  very  approximate.  An 


Kinetic  1 

Table  4.2 
energy  e-Polding 

Dissipation 

Time  Scale, 

for  Nixed  Layer 

Depth  d  >  20  m 

N0k  (s’*1nr1, 

Case  T( inertial  periods) 

3.125  x  1(H 

5 

30 

6.25  x  10“7 

6 

8 

1.25  x  10“* 

1 

2.4 

2.50  x  10“* 

7 

1 

approximate  empirical  formula  for  the  dissipation  time  scale 
is 

T( inertial  periods)  -  (2.175  x  10“6/NQk)1,7  (4.11) 

This  formula  is  correct  only  for  the  particular 
value  used  in  Cases  1,  5,  6,  and  7  for  dimensionless  mixed 
layer  depth  d  ■  20  m.  Since  internal  waves  do  not  feel  the 
computational  bottom,  the  only  physically  meaningful  vertical 
length  scale  is  the  mixed  layer  depth.  Therefore,  the 


Depth-Averaged  Kinetic  Energy  Density  (kg 


Case  5 


Time  (inertial  periods) 

Figure  4.6  Depth-averaged  kinetic  energy  in  the  surface  layer 
for  Cases  5,  6,1,  7,  with  NqK  *  3.125  x  10“7, 

6.25  x  10“7,  1.25  x  10-6,  and  2.5  x  10"6  s'^-tir1, 
respectively. 


relevant  dimensionless  parameter  is  Nkd/f,  where  d  is  the 
mixed  layer  depth.  Using  this  parameter,  we  can  generalize 
(4.11)  to  obtain 


T(days ) 


1  T  sine  Mm)  1  1.7 
TTn"§  [260  N  (cph)d(m)  J  ' 


(4.13) 


where  6  is>  the  latitude  and  X  »  2 n/k  is  the  horizontal 

wavelength.  Table  4.3  shows  a  tabulation  of  (4.13),  for 
representative  values  of  X  and  N0.  Blank  areas  in  the 
table  are  outside  the  range  covered  by  Cases  1,  5,  6,  and  7. 


Table  4.3 

Kinetic  Energy  e-Folding  Dissipation  Time  Scale  (days). 


Latitude 

e  «  40 *N,  Mixed 

Layer 

Depth  d 

■  20  m 

X(  km) 

25 

50 

100 

200 

2 

3.3 

10,6 

N  (cph) 

4 

1 .0 

3.3 

10.6 

o 

8 

0.3 

1 .0 

3.3 

10.6 

16 

0.3 

1 .0 

3.3 

The  scale  X  is  the  wavelength  associated  with  an 
atmospheric  front.  In  a  case  study  of  a  cold  front  by  HOBBS 
et  al.  (1980),  the  wavelength  of  the  surface  wind  veering  is 
about  120  km.  Squall  lines  and  narrow  cold  fronts  are 
observed  to  have  associated  wavelengths  as  short  as  10  km 
(CARBONE,  1982;  HOBBS  and  PERSSON,  1982).  Direct  estimates 
from  observations  of  inertial  oscillations  show  X  to  vary 
over  a  wide  range.  THOMSON  and  HUGGETT  (1981)  estimated  X  in 


Section  5 

NONUNIFORM  STRATIFICATION 


Z*  *  d  -  Z. 


(5.1) 


Figure  5.2  shows  the  response  of  the  v-component  of 
velocity  for  Case  1,  where  no  bottom  porosity  damping  is 
included.  Thus,  reflections  off  the  bottom  of  the  computa¬ 
tional  domain  are  allowed.  The  z-axis  is  unstretched,  and 
ranges  from  0  to  D.  The  time  axis  ranges  from  t  ■  0  to 


We  can  extend  the  results  of  Section  4  by  consid¬ 
ering  a  somewhat  more  realistic,  nonuniform  distribution  for 
N(z); 


z’  <  0, 


N(z')  */ Nj  (e/c)  z'exp(-z'/c) 

^(N^/a)  [ ( a—  1 )  (e/c)  z'exp(-z'/c)  +  l] 


0  <  z'  <  c. 


c  <  z'  , 


This  profile  is  shown  in  Figure  5.1.  N(z)  is  zero  in  the 
surface  layer  (of  depth  d) ,  has  a  peak  value  Nx  at  a 
depth  c  below  the  base  of  the  surface  layer,  and  decays 
asymptotically  to  N^/a  with  large  depth.  The  eddy 
diffusivity  profile  is  the  step  function  given  by  (4.2),  and 
is  also  shown  in  Figure  5.1.  As  in  Section  4,  we  again  apply 
a  surface  stress  impulse  at  time  t  *  0  in  the  x-direction. 
Table  5.1  lists  the  parameters  for  all  cases  discussed  in 
this  section. 

5.1  SOLUTION  OVERVIEW 


t 


Case 

1  6.25  x 

2  6.25  x 

3  6.25  x 

4  1.25  x 

5  2.5 

6  5 

7  5 

8  5 

9  5  x  10~6  10 


Porosity 

No 

Yes 

Yes 

Yes 

Yes 

Yes 


2.4 

0.02 

Yes 

4.8 

0.02 

Yes 

9.6 

0.02 

Yes 

Table  5.1 

Parameters  for  Cases  in  Section  5 


N-jk  (s^m"1) 

d(m) 

c(m) 

(m2s_1) 

6.25  x  10~7 

80 

40 

0.16 

6.25  x  10“7 

80 

40 

0.16 

6.25  x  10“7 

20 

10 

0.02 

1.25  x  10“6 

20 

10 

0.02 

2.5  x  10"6 

20 

10 

0.02 

VO 

1 

O 

X 

in 

20 

10 

0.02 

x  10~6  10 

x  Ifl'6  10 


r 


I 
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line  axis  runs  from  0  to  4  inertial  periods 


8n/f,  or  4  inertial  periods.  Momentum  becomes  evenly  dif¬ 
fused  through  the  surface  layer  during  the  first  1/4  inertial 
period.  The  flow  in  the  surface  layer  oscillates  at  a  near- 
inertial  frequency.  Internal  waves  propagate  into  the 
stratified  region,  and  a  velocity  maximum  develops  just  below 
the  surface  lay«r.  These  features  are  qualitatively  similar 
to  the  uniformly  stratified  case  shown  in  Figure  4.1,  but 
since  the  stratification  in  Figure  5.2  is  weaker,  the 
vertical  phase  velocity  is  greater  and  the  bands  of  constant 
phase  are  more  nearly  vertical.  After  time  t  =  4n/f  (2 
inertial  periods),  reflections  off  the  bottom  boundary  become 
important  and  generate  an  interference  pattern.  We  show  this 
solution  in  order  that  we  can  make  a  comparison  with 
observations  in  a  shallow-bottom  sea,  in  Section  5.3. 

In  the  deep  ocean,  the  ocean  bottom  has  a 
negligible  effect  on  the  propagation  of  near-inertial  waves. 
Figure  5.3  shows  the  v-component  of  velocity  for  Case  2, 
where  a  bottom  porosity  distribution  has  been  introduced  in 
the  approximate  range  0  £  z  <  0.35D.  The  z-axis  here  ranges 
from  0.35D  to  D,  (the  depth  ranges  from  the  surface  down  to 
520  m)  in  order  to  avoid  overlapping  this  porosity  distribu¬ 
tion.  Internal  waves  are  absorbed,  and  thus  the  model 
simulates  an  ocean  of  infinite  depth.  The  principal  dif¬ 
ference  between  Figures  5.2  and  5.3  is  that  Figure  5.3  lacks 
the  interference  pattern  in  the  deeper  half  of  the  computa¬ 
tional  domain,  which  results  from  bottom  reflections. 

5.2  SURFACE  LAYER  DISSIPATION 

Close  inspection  of  Figures  5.2  and  5.3  shows  an 
interesting  feature  which  is  absent  from  the  uniformly 
stratified  solution.  The  amplitude  of  the  surface  layer 
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inertial  oscillations,  as  well  as  the  amplitude  of  the 
velocity  maximum  just  below  the  surface  layer  are  both 
periodic,  and  about  180*  out  of  phase  from  each  other.  This 
behavior  is  analogous  to  that  of  two  coupled  oscillators. 
This  "frequency  interference  effect"  results  from  the  peak  in 
the  Vaisala  frequency  profile,  which  causes  certain  internal 
wave  components  to  resonate  and,  in  a  sense,  to  become 
partially  trapped.  Figure  5.4  shows  the  oscillations  of 
depth-averaged  kinetic  energy  density  in  the  surface  layer 
0.9  £  z/D  £  1.0  and  in  the  strongly  stratified  layer  0.8  £ 
z/D  £  0.9,  for  Case  2.  The  period  of  these  oscillations  is 
about  4.6  inertial  periods,  and  remains  constant  throughout 
the  integration. 

Figure  5.5  shows  the  time  evolution  of  depth 
averaged  kinetic  energy  density  in  the  surface  layer  for 
Cases  3-6.  The  dissipation  rate  increases  with  increasing 
Ni .  The  e-folding  dissipation  time  scale  is  tabulated  in 
Table  5.2.  These  time  scale  estimates  are  very  approximate. 
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Table  5.2 

Kinetic  Energy  e-Folding 

Dissipation 

Time  Scale,  for 

Nixed  Layer 

Depth  d  «  20  m 

Njk  (s”1!!!”1) 

Case  T( inertial  periods) 

6.25  x  10"* 

3 

13 

1.25  x  10“6 

4 

3.5 

2.5  x  10"® 

5 

1.25 

5  x  10"6 

6 

0.42 

Depth- Averaged  Kinetic  Energy  Density  (kg 


Time  (Inertial  Periods) 


Fig.  5.4  Depth  averaged  kinetic  energy  in  the  surface  layer 
0.9  z/d  £  1.0,  and  in  the  stratified  interior 
layer  0.8  <  z/D  <  0.9,  for  Case  2. 


Depth-Averaged  Kinetic  Energy  Density  (kg 


An  approximate  empirical  formula  for  the  dissipation  time 
scale  is 


T(inertial  periods)  ■  ( 3  x  10”6/N1 k) 1 *7.  (5.2) 

The  dissipation  time  scale  given  by  (5.2)  is  greater  than 
that  in  (4.11),  but  we  remind  the  reader  that  N0  in  (4.11) 
refers  to  a  uniform  stratification,  while  in  (5.2) 

refers  to  the  peak  level  of  stratification. 

As  in  Section  4.4,  we  can  express  a  generalized 
formula  for  the  dissipation  time  scale,  in  units  of  days,  as 
follows: 

T(days)  *  sfs*  (Cp£">dUj]  *  (5*3) 

Table  5.3  is  a  tabulation  of  (5.3),  for  represen¬ 
tative  values  of  \  and  Nj .  Blank  areas  in  the  table  are 

out  of  the  parameter  range  covered  by  Cases  3-6.  Typical 

values  for  Nj  are  of  order  20  cph.  with  this  value, 
inertial  oscillations  produced  by  wind  fields  with  dominant 
wavelengths  less  than  200  km  are  all  dispersed  from  the 
surface  layer  within  a  few  days. 

These  results  cannot  be  readily  generalized  to 

Vaisala  frequency  profiles  with  parameters  different  from 
those  listed  for  Cases  3-6  in  Table  5.1.  The  strength  of  the 
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"frequency  interference  effect"  depends  on  the  relative 
amplitude  of  the  peak  of  the  vaisala  frequency  distribution 


Table 

5.3 

Kinetic  Energy 

e-Folding  Dissipation  Time  Scale 

(days) , 

for  Latitude 

6  ■  40  N  and 

Mixed  Layer 

Depth  d 

■  20m 

k  (km) 

25 

50 

100 

200 

5  1.0 

3.3 

10.6 

N1  (cph) 

10  0.3 

1.0 

3.3 

10.6 

20 

0.3 

1.0 

3.3 

30 

0.5 

1.6 

(the  parameter  a  in  (5.1)).  When  a»1,  the  problem  reverts 
back  to  the  uniform  stratification  profile  discussed  in 
Section  4.  As  a  increases ,  the  "frequency  interference 
effect"  becomes  progressively  stronger. 

The  resonant  frequency  of  the  interference  effect 
depends  on  the  magnitude  N;  k,  as  well  as  on  the  pycnocline 
thickness  parameter  c.  The  dependence  of  the  resonant 
frequency  «*>r  on  Nik  can  be  estimated  for  the  particular 
cases  illustrated  in  Figure  5.5 

«r  »  f [  1  +  (2.4  x  105  N,k) 2]  ^2.  (5.4) 

The  dependence  of  the  resonant  frequency  on  the  pycnocline 
thickness  is  presented  in  Figure  5.6,  which  shows  the  time 
evolution  of  depth  averaged  kinetic  energy  in  the  surface 
layer  for  Cases  7-9.  These  cases  are  identical  except  in 
their  pycnocline  thickness  parameter  c.  The  dissipation  rate 


Time  (inertial  periods) 


Figure  5 . 


6  Depth-averaged  kinetic  energy  density 
in  the  surface  layer  for  Cases  7-9 , 
with  pycnocline  thickness  parameter 
c  »  2.4,  4.8,  9.6  m. 


initially  increases  with  increasing  c.  The  reason  is  that 
with  c  increasing,  the  high  stratification  region  extends  to 
a  greater  depth,  and  allows  faster  vertical  dispersion. 
After  some  time  has  elapsed,  the  wave  component  whose 
vertical  wavelength  is  resonant  with  the  peak  in  the  Vaisala 
frequency  distribution  reaches  the  strong  pycnocline,  and 
interferes  with  the  dominant  inertial  oscillations  in  the 
surface  layer.  The  frequency  of  this  interference  effect 
increases  with  increasing  pycnocline  thickness.  The 
amplitude  of  the  interference  effect  (that  is,  the  amplitude 
of  the  kinetic  energy  oscillation  in  absolute  terms)  is 
independent  of  the  pycnocline  thickness,  c.  After  the 
resonant  wave  component  has  had  time  to  propagate  into  the 
uniformly  stratified  region,  the  dissipation  of  surface  layer 
kinetic  energy  continues,  at  a  rate  independent  of  c. 

5.3  MODEL-DATA  COMPARISON:  BALTIC  SEA 

In  this  section  we  present  a  qualitative  comparison 
between  the  model  results  and  current  meter  observations  in 
the  Baltic  Sea.  The  observations  were  taken  at  56*05. 7'N, 
18*44.4' E,  a  location  where  the  bottom  depth  is  105  m,  and 
are  described  by  KRAUSS  (1981).  Current  meters  were  deployed 
at  10  m  intervals  between  the  depths  30  and  100  m. 

Figure  5.7  shows  the  u-component  of  velocity  from  8 
September  to  21  September  1977,  from  KRAUSS  (1981).  Three 
storms  occurred  during  this  interval,  a  relatively  weak  one 
with  wind  speeds  up  to  22  m/sec  on  10  September,  and  two 
stronger  storms  with  wind  speeds  up  to  29  m/sec  on  12  and  14 
September.  We  wish  to  concentrate  our  attention  on  the 
interval  from  12  to  16  September.  During  this  interval,  the 
base  of  the  mixed  layer  deepened  from  about  25  to  30  m,  as 
shown  in  the  smoothed  temperature  distribution  in  Figure  5.8. 
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Figure  5.7  u-component  of  velocity  in  the  Baltic  Sea 
during  the  time  interval  8-20  September, 

1977,  from  KRAUSS  (1981).  Contour  intervals 
are  plotted  for  every  5  cm/sec.  Heavy  line 
is  zero  velocity;  full  lines  are  positive  and 
dashed  lines,  negative. 


The  storm  on  12  September  increased  the  amplitudes 
of  the  inertial  waves,  while  the  phase  of  the  storm  on  14 
September  was  such  as  to  reduce  the  amplitudes  somewhat. 
During  the  interval  12  to  16  September  there  is  a  definite 
trend  toward  smaller  vertical  phase  velocities.  In  other 
words,  the  bands  of  constant  phase  increasingly  slope  away 
from  the  vertical.  This  feature,  plus  the  abrupt  change  in 
phase  from  the  mixed  layer  to  the  thermocline  (seen  more 
clearly  in  Krauss'  Figure  11,  but  not  shown  here)  are 
strikingly  similar  to  the  model  results  in  Figure  5.2. 

In  a  shallow  bottom  sea,  bottom  reflections  become 
important  several  inertial  periods  after  a  wind  event  has 
occurred.  From  about  15  to  16  September  (as  well  as  on  19 
September)  there  are  indications  of  bottom  reflections. 
These  reflections  show  up  as  a  sharp  phase  shift  at  about 
70  m  depth.  This  phase  shift  is  similar  to  that  which 
develops  in  Figure  5.2  (Case  1)  about  3  inertial  periods 
after  the  wind  impulse. 


5.4  MODEL-DATA  COMPARISON:  DEEP  OCEAN 

In  this  section  we  present  a  qualitative  comparison 
between  our  model  results  and  current  observations  in  the 
deep  ocean.  For  this  comparison,  we  consider  the  Mixed  Layer 
Experiment  (MILE)  which  was  carried  out  in  the  northeastern 
Pacific  Ocean  near  Ocean  Weather  Station  P.  The  MILE-1 
mooring,  located  at  49*37'N,  145*6'W  was  connected  to  a 
surface  float  in  4360  m  of  water.  A  complete  description  of 
the  MILE  setting  may  be  found  in  DAVIS  et  al.  (1981).  The 
MILE-1  mooring  collected  data  over  a  20-day  period,  from  19 
August  to  7  September.  Current  meters  collected  data  samples 
at  a  rate  of  32  hr”1,  which  we  filtered  to  1  hr"1. 


Figure  5.9  shows  a  contour  plot  of  the  u  component 
of  velocity.  As  the  two  uppermost  current  meters  delivered 
incomplete  records,  the  depth  range  covered  by  these  plots  is 
from  11  to  121  m.  Inertial  oscillations  dominate  much  of  the 
20-day  records.  The  inertial  period  at  this  latitude  is 
about  17  hr.  Vertical  shear,  marked  by  an  abrupt  change  in 
phase,  is  concentrated  near  the  top  of  the  thermocline,  at  a 
depth  of  about  40  m.  Below  40  m,  the  bands  of  constant  phase 
are  for  the  most  part  nearly  vertical.  The  most  significant 
exception  is  from  29-30  August,  when  the  phase  progression  is 
downward.  The  implication  is  that  the  propagation  direction 
of  a  local  internal  wave  packet  is  upward. 

Figure  5.10  shows  the  wind  stress  components  during 
MILE ,  reproduced  from  DAVIS  et  al.  (1981).  Two  major  wind 
events  occurred,  on  21-24  August  and  30  August- 1  September. 
Both  of  these  wind  events  generated  inertial  oscillations  in 
the  mixed  layer,  during  22-23  August  and  1-3  September, 
respectively. 

Both  in  the  mixed  layer  and  below,  the  inertial 
oscillations  are  intermittent,  and  decay  within  about  2  days 
after  their  generation.  Their  peak  amplitudes  are  about 
20-30  cm  sec'1.  The  strongest  peak  amplitudes  in  the  mixed 
layer  occur  on  20  August,  23  August,  25  August,  2  September, 
and  5  September,  and  the  strongest  peak  amplitudes  in  the 
thermocline  occur  on  19  August,  21-22  August,  24  August,  and 
6  September.  There  seems  to  be  a  2-3  day  cycle,  in  which 
strong  inertial  waves  oscillate  between  the  mixed  layer  and 
upper  seasonal  thermocline.  This  oscillation  is  similar  in 
character  to  the  model-predicted  frequency  interference 
effect,  described  in  Section  5.2.  The  peak  amplitudes 
which  occur  in  the  mixed  layer  on  20  August,  25  August,  and 
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5.10  Wind  stress  to  the  north  (upper 
trace)  and  to  the  east  (lower 
trace)  at  Station  P  from  19 
August  to  7  September,  1977. 
From  DAVIS  et  al.  (1981). 
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5  September  do  not  seem  to  result  directly  from  wind  events. 
This  observation  lends  support  to  the  hypothesis  that  the 
frequency  interference  effect  is  responsible  for  these 
oscillations. 


Section  6 

SUMMARY  AND  RECOMMENDATIONS 


During  FY81,  we  developed  a  linear  model  of  iner¬ 
tial  waves  in  a  shallow  sea  (Rubenstein  1981b,  1982).  This 
model  predicted  the  near-inertial  frequency  response  of  the 
ocean  to  wind  events.  The  model  equations  were  solved  with  a 
Fourier-Chebyshev  expansion  in  the  spatial  dimensions,  and 
numerically  in  time.  Delta-function  wind  stress  impulses 
were  applied  at  the  surface.  The  response  in  the  shallow  sea 
was  shown  to  be  governed  by  vertical  modes.  The  frequencies 
of  these  modes  satisfied  the  internal  wave  dispersion  rela¬ 
tion. 


During  FY82,  we  extended  this  model  to  the  deep 
ocean,  as  described  in  this  report.  While  the  equations  are 
very  similar  to  those  of  KRAUSS  (1978,  1979),  we  have  intro¬ 
duced  a  porous  damping  region  which  minimizes  reflections 
from  the  bottom  of  the  computational  domain.  Thus,  the  model 
simulates  an  ocean  of  "infinite"  depth.  Vertical  modes 
associated  with  standing  waves,  prevalent  in  our  earlier 
solutions,  do  not  develop  in  this  model.  We  successfully 
validated  the  model  against  analytic  solutions  of  the 
diffusion  problem. 

We  applied  this  model  to  determine  the  rate  at 
which  inertial  oscillations  in  the  surface  layer  are  dissi¬ 
pated.  We  looked  first  at  a  uniformly  stratified  ocean. 
High  frequency  internal  wave  components  had  the  fastest  group 
velocities,  and  immediately  after  a  wind  event  they  left  the 
surface  layer  soonest.  Lower  frequency,  near-inertial  waves 
were  left  behind,  and  propagated  downward  much  more  slowly. 
The  dissipation  rate  of  surface  layer  kinetic  energy  is  a 
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function  of  the  ratio  Nd/fX,  where  X  is  the  horizontal  wave¬ 
length  of  the  wind  event  (a  front  or  storm  system),  d  is  the 
mixed  layer  depth,  f  is  the  local  inertial  frequency,  and  N 
is  the  dimensional  Vaisala  frequency,  a  measure  of  the 
stratification. 

We  also  considered  a  more  realistic  nonuniform 
stratification  distribution,  with  a  strong  inviscid  pycno- 
cline  immediately  below  a  diffusive  surface  layer.  For  cases 
where  the  peak  Vaisala  frequency  was  of  order  20  cph,  and 
the  wavelength  associated  with  the  wind  event  was  less 
than  200  km,  the  kinetic  energy  in  the  surface  layer  dissi¬ 
pated  within  a  few  days.  This  result  is  in  agreement  with 
simple  mixed  layer  inertial  oscillation  prediction  models 
(POLLARD  and  MILLARD,  1970;  KUNDU ,  1976;  POLLARD,  1980). 
These  models  show  the  best  agreement  with  observational  data 
when  their  damping  factors  correspond  to  e-folding  decay 
times  (for  the  amplitude  of  inertial  oscillations)  in  the 
range  2  to  8  days.  Therefore,  vertical  dispersion  alone  is 
adequate  to  account  for  the  observed  dissipation  rates. 

When  a  sufficiently  peaked  pycnocline  is  present,  a 
frequency  interference  effect  occurs.  A  particular  low  fre¬ 
quency  wave  component  resonates  within  the  pycnocline,  and 
interferes  with  the  dominant  inertial  oscillations  in  the 
surface  layer.  The  strength  of  the  interference  depends  on 
the  peakedness  of  the  pycnocline.  The  resonant  frequency 
depends  on  the  stratification  and  thickness  of  the  pycno¬ 
cline.  As  a  result  of  this  effect,  the  kinetic  energy  in  the 
surface  layer  may  not  dissipate  monotonically.  Energy  may  be 
alternately  exchanged  between  the  surface  layer  and  the 
pycnocline  immediately  below.  We  examined  deep-ocean 
observations  taken  during  MILE,  and  found  some  evidence  of 


this  frequency  interference  effect.  Several  of  the  incid¬ 
ences  of  inertial  oscillations  in  the  mixed  layer  did  not 
seem  to  result  directly  from  any  particular  local  wind  event. 
Instead,  kinetic  energy  was  cyclically  exchanged  between  the 
surface  layer  and  the  upper  seasonal  thermocline.  The  cycle 
period  was  2  to  3  days. 

We  also  compared  model  simulations  with  current 
meter  observations  in  the  shallow  (105  m)  Baltic  Sea.  We 
found  that  vertical  phase  velocity  generally  decreases  after 
the  occurrence  of  a  storm.  Our  comparisons  reveal  that 
bottom  reflections  become  important  after  a  few  inertial 
periods  have  elapsed  after  a  storm. 

A  number  of  ocean  mixed -layer  prediction  models 
incorporate  damping  factors  which  parameterize  the  dispersion 
of  near-inertial  internal  waves  (POLLARD  and  HILLARD/  1970; 
KUNDU,  1976;  CLANCY  et .  al . ,  1981).  The  magnitudes  of  these 
damping  factors  are  known  to  vary  with  changing  environmental 
conditions,  but  are  determined  through  the  method  of  obtain¬ 
ing  a  "best  fit"  with  a  data  set.  The  results  of  the  present 
study  should  allow  investigators  to  parameterize  the 
dispersion  of  internal  waves  with  better  accuracy. 

One  such  mixed-layer  model,  known  as  TOPS  (Thermo¬ 
dynamical  Ocean  Prediction  System)  was  developed  at  the  Naval 
Ocean  Research  and  Development  Activity  for  operational  use 
at  the  Navy's  Fleet  Numerical  Oceanography  Center.  TOPS  is 
used  to  predict  the  ocean  mixed-layer  environment  over  large 
geographical  regions.  TOPS  incorporates  damping  factors 
which  determine  the  rate  of  attenuation  of  inertial  oscilla¬ 
tions  in  the  mixed-layer.  Using  the  results  of  the  present 


study,  it  may  be  possible  to  construct  a  more  realistic 
parameterization  of  the  attenuation  rate,  which  depends  on 
the  local  environment. 

The  model  described  in  this  report  is  linear;  as  a 
result,  the  profiles  of  eddy  diffusivity  and  Vaisala  fre¬ 
quency  are  constant  in  time.  It  would  be  more  realistic  to 
include  turbulent  buoyancy  fluxes  and  to  allow  the  Vaisala 
frequency  profile  to  evolve  with  time.  The  eddy  diffusivity 
could  be  parameterized  in  terms  of  the  local  shear,  or 
Richardson  number.  We  recommend  that  the  current  model  be 
extended  to  include  these  changes.  In  this  way,  we  should 
obtain  better  insight  into  the  dynamics  of  the  upper  ocean  in 
the  aftermath  of  the  passage  of  wind  events  and  storms. 
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